1 Introduction

Most of Central America—here loosely defined as the continental land between southern México and Panamá—experienced below average rainfall between 2015 and 2019 (Fig. 1a, Fig. S1). This prolonged rainfall deficit resulted in severe and continued drought conditions affecting the livelihoods of millions of people living in Guatemala, El Salvador, Honduras, Nicaragua, and Costa Rica (Palencia 2014; Moloney 2018; Masters 2019). Specifically, years 2019, 2015, and 2018 have seen the largest (in magnitude) negative anomalies of summer precipitation in Central America over the last century (Fig. 1b). Furthermore, the period 2015–2019 experienced the lowest 5-year mean May-September (MJJAS) precipitation (Fig. 1c) over the last century.

Fig. 1
figure 1

Recent multi-year meteorological drought in Central America. a Map of 2015–2019 May-to-September precipitation anomaly (GPCC dataset; Schneider et al. 2013) over Central America, relative to the 1921–1970 climatological mean. The grey outline denotes the study region used to define the Summer Rainfall Index (SRI). b Histogram of annual MJJAS SRI anomalies and c 5-yr mean MJJAS SRI anomalies. d Five-year running mean of MJJAS SRI anomalies for four observational products (GPCC, CRU, CMAP, and CHIRPS, see Section 2.1 for details on these datasets). The horizontal dashed line denotes the lowest observed 5-year anomaly (2015–2019)

Full size image

The 2015–2019 exceptionally dry period is at the end of an approximately 40-year time period (from the late 70s’ to present) characterized by a negative trend in summertime precipitation (Fig. 1d, Fig. 2, see also, e.g., Neelin et al.2006). Summer precipitation over Central America shows a marked decadal variability resulting from modes of large scale variability like El Niño-Southern Oscillation, the Atlantic Multidecadal Oscillation, and the North Atlantic Oscillation (Giannini et al. 2000, 2001; Wang 2007; Hastenrath and Polzin 2013; Muñoz-Jimènez et al. 2019; Anderson et al. 2019; Hidalgo et al. 2019). Climate projections from global and regional models consistently suggest a reduction of summertime precipitation between May and September in Central America by the end of this century, along with higher summer temperatures (Giorgi 2006; Rauscher et al. 2008; Fuentes-Franco et al. 2015; Hidalgo et al. 2013; Almazroui et al. 2021; Depsky and Pons 2020). This raises the question of whether Anthropogenic Climate Change (ACC) has contributed to the 2015–2019 meteorological drought and to the 40-year negative trend, increasing their likelihood.

Most of Central America features a monsoonal climate, with the summertime wet season characterized by a bimodal regime with two precipitation peaks in June and September separated by the midsummer drought (MSD; Magaña et al. 1999; Gamble et al. 2008, see also Fig. S1). In areas of Central America on the side of the Caribbean Sea, the precipitation distribution is more uniform during the wet season though, with higher annual accumulations compared to the rest of Central America (Alfaro 2002; Martinez et al. 2019). Because of the highly seasonal character of Central American rainfall, even small changes in rainfall can have a substantial impact on regional water resources and local rain-fed agriculture. Indeed, crop shortages due to extreme drought conditions over the last decade and the subsequent impoverishment have been identified as the main driver of mass migrations from Central America in recent years (World Food Programme 2020). Climate migration from Central America is expected to continue and ramp up in the next decades due to stronger impacts of global warming combined with rapid population growth (Rigaud et al. 2018; Hidalgo and Alfaro 2012).

Observational evidence of more prolonged MSDs along with more dry days has been recently discussed in the literature (Maurer et al. 2017; Anderson et al. 2019), but there has been little progress so far in attributing these changes to ACC. Herrera et al. (2018) examined the 2013–2016 Pan-Caribbean drought, which affected mostly the Caribbean islands, and concluded that ACC played a major role in increasing the drought’s severity through the effects of higher temperature on potential evapotranspiration. This made its attribution to ACC somewhat less problematic, given the observed, long-term warming of the region (Aguilar et al. 2005; Pachauri and Meyer 2014). Attribution of a meteorological drought to ACC is instead more troublesome, since this is determined by local and large-scale circulation changes which, in turn, are generally less directly controlled by temperature (NAS 2016). We have no evidence of specific attribution studies connecting ACC to the recent 2015–2019 meteorological drought in Central America. Neelin et al. (2006) analyzed the recent multi-decadal summer drying trend over Central America and concluded that the observed trend cannot be unambiguously attributed to ACC, that is, it cannot be excluded that it is due to natural interdecadal variability.

This study has two objectives. First, to determine if the recent drying trend observed from the late 1970s is attributable to ACC. Second, to evaluate whether the 2015–2019 Central American meteorological drought has been influenced, i.e., made more likely, by ACC. These two questions will be addressed from a meteorological point of view by looking at the impact of climate change on the precipitation deficit. We realize that this is only the first step in analyzing the drivers of a complex phenomenon like a drought, which is not only caused by large scale circulation anomalies but also by complex interactions of local precipitation, wind, temperature, land use, and water management (Mishra and Singh 2010; AghaKouchak et al. 2021).

We will use five different large ensembles—including two high-resolution (i.e., 0.5∘ horizontally) ensembles—to estimate the contribution of ACC and distinguish it from that of internal variability. Large ensembles of climate models (Deser et al. 2020) are very useful for climate risk and attribution studies (Otto et al. 2018; Swain et al. 2018; Pascale et al. 2020) as they provide thousands of years of data and thus allow for a direct reconstruction of the underlying probability distribution of hydroclimatic extremes without relying on a hypothesized statistical model of extremes (Van der Wiel et al. 2019). Central American summer rainfall has large magnitude interannual and decadal variability (Giannini et al. 2000, 2001; Wang 2007; Hastenrath and Polzin 2013; Muñoz-Jimènez et al. 2019; Anderson et al. 2019) and a large ensemble is, thus, also a powerful method to isolate, at the decadal timescale, internal variability from the forced signal (Kay et al. 2015; Deser et al. 2020).

2 Data and methods

2.1 Observation data

To take into account uncertainty in observed precipitation over Central America, four different observational precipitation datasets are used in this study: (1) the Global Precipitation Climatology Centre (GPCC) dataset version 7, at 0.5∘ horizontal resolution (Schneider et al. 2013), (2) the Climate Research Unit high-horizontal resolution grids of monthly rainfall at the University of East Anglia, version 3.24, at 0.5∘ horizontal resolution (Harris et al. 2013), (3) the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS), at 0.05∘ horizontal resolution (Funk et al. 2015), and (4) the CPC Merged Analysis of Precipitation (CMAP), at 2.5∘ (Xie and Arkin 1997) horizontal resolution.

For atmospheric variables, like, e.g., the 925 hPa winds, we employ the ERA5 reanalysis (Hersbach et al. 2019). ERA5 is the most recent reanalysis product from the European Centre for Medium-Range Weather Forecast and it features major improvements such as higher horizontal and vertical resolution (0.28∘ and 139 vertical levels), an updated and enhanced version of the Integrated Forecast System (IFS) Earth system model and associated observational assimilation system (Hersbach and Dee 2016).

Sea surface temperatures (SSTs) are taken from the Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5) dataset (Huang et al. 2017), which is a global monthly SST dataset from 1854 to present on a 2∘× 2∘ grid derived from the International Comprehensive Ocean-Atmosphere Dataset.

2.2 Modelling data

We use a suite of large ensemble simulations from the the S eamless System for P rediction and EA rth System R esearch (SPEAR; Delworth et al., 2020). SPEAR is the newest modeling system for seasonal to multi-decadal prediction developed at NOAA Geophysical Fluid Dynamics Laboratory (GFDL) and shares underlying component models with the CM4 climate model (Held et al. 2020). The SPEAR atmospheric model is run at different horizontal resolutions (atmosphere, land) in this paper: 0.5∘(SPEAR_MED) and 1∘ (SPEAR_LO), and it has 33 atmospheric levels in the vertical. More details about SPEAR and the SPEAR large ensemble can be found in Delworth et al. (2020), Lu et al. (2020), Pascale et al. (2020), and Murakami et al. (2020). The SPEAR_MED large ensemble is characterized by a horizontal grid-spacing that is finer than those of most other available large ensembles, which makes the SPEAR_MED ensemble an unprecedented and unique tool to study regional climates.

To evaluate forced vs. natural variability, we use four different numerical experiments: (1) CTRL, a long-term control run with constant preindustrial (1850) forcing; (2) NATURAL, an ensemble driven only by natural forcing (i.e., volcanic eruptions and solar cycles); (3) ALLFORC4.5, an ensemble driven by observed anthropogenic and natural forcing up to 2014 (HIST), and then according to the Shared Socioeconomic Pathway (SSP2-4.5) developed for the Coupled Model Intercomparison Project Phase 6 (CMIP6) (O’Neill et al. 2017; Eyring et al. 2016; Riahi et al. 2017); and (4) ALLFORC8.5, an ensemble driven by observed anthropogenic and natural forcing up to 2014 (HIST) and then according to the Shared Socioeconomic Pathway (SSP5-8.5). The SSP5-8.5 (“high emissions”) represents the high end of the range of future scenarios (radiative forcing by 2100 of 8.5 W m2) and it updates CMIP5 RCP8.5 (Kriegler et al. 2017). SSP5-8.5 is compatible with a future development heavily based on fossil fuels. The SSP2-4.5 (“middle of the road emissions”) instead sits in the middle of the range of future forcing pathways and updates the CMP5 RCP4.5 (O’Neill et al. 2017).

We also analyze four additional large ensembles to assess model uncertainty: SPEAR_LO, the Forecast-Oriented Low Ocean Resolution model with flux adjustment, FLOR_FA (Zhang and Delworth 2018), the Community EARTH System Model Large Ensemble, CESM-LENS (Kay et al. 2015), and the Max Planck Institute Grand Ensemble, MPI-GE (Maher et al. 2019). While SPEAR_LO forced runs follow the CMIP6 SSP5-8.5, the last three ensembles are available with various CMIP5 scenarios. Table 1 provides additional details of the SPEAR and other models’ experiments used in this study. Although developed following different protocols (CMIP5 vs CMIP6) and thus being different in the details of the future anthropogenic forcings, SSP5-8.5 and SSP2-4.5 are comparable to RCP8.5 and RCP4.5, respectively, when evaluated on the base of carbon dioxide emissions and radiative forcing (Kriegler et al. 2017; Meinshausen et al. 2500).

Table 1 Description of models and simulations used in this study
Full size table

2.3 Model evaluation

An important step in multi-model analysis of regional hydroclimates is to evaluate each model’s performance to determine which models more reliably represent the hydroclimate of a certain region. We assessed each model’s performance in terms of (1) the annual cycle of precipitation over Central America (Fig. S3-S4 and Table 2), (2) the amplitude of the interannual, multi-annual, and decadal variability of the Central American summer rainfall (Fig. S5) and ability to simulate the tail of the probability distribution of the 5-year MJJAS rainfall anomalies (that is, multi-annual negative rainfall deficits leading to droughts like the 2015-2019 one, Fig. S6), and (3) the ability to capture the remote SST drivers of Central American rainfall in the Atlantic and Pacific Ocean (Fig. S7).

Table 2 Evaluation of the large ensemble models
Full size table

SPEAR_MED has the most realistic representation of the Central American summer rainfall patterns (Fig. S3) and seasonal cycle (Fig. S4) according to six different precipitation metrics (Table 2), followed by SPEAR_LO and FLOR_FA. All models reproduce quite well the magnitude of rainfall variability, with the exception of CESM that underestimates it (Fig. S5). The SPEAR models best reproduce the left tail of the 5-year MJJAS rainfall anomaly cumulative probability distribution (Fig. S6). All models capture—though exaggerating it—the remote influence of the Atlantic and Pacific oceans on Central American rainfall (Fig. S7), with MPI-GE and FLOR_FA exaggerating the influence of the tropical Pacific and partially missing that of the eastern tropical Atlantic. Collectively, these results suggest that all models used in this study are appropriate tools for use in characterizing future changes in regional precipitation over Central America, with the SPEAR_MED and SPEAR_LO the best performing models.

Finally, we also evaluate if the models’ historical trends of MJJAS rainfall are consistent with observations over Central America over the period 1979–2019 (Fig. S8). To do so, we compute rainfall trends over the period 1979–2019 in GPCC and compared them with individual members of the ALLFORC8.5 ensembles over the same time period. If the observed trend at one grid point was within the range of those simulated by the 30 ensemble members, then we said that the model is consistent with observations in that grid box. We found that models are consistent with observations over most of the grid points (Fig. S8).

2.4 Climatological indexes and event definition

As precipitation in Central America is mostly concentrated from May to September, we use a Central American Summer Rainfall Index (SRI hereafter) to monitor the interannual and monthly rainfall anomalies. The SRI is defined as the area-averaged MJJAS precipitation over all land grid points between 11∘N–18∘N and 95∘W–84∘W (Fig. 1a). This is the area most affected by the 2015–2019 meteorological drought (Fig. 1a), and will therefore be our study region.

A more intense Caribbean Low-Level Jet (CLLJ) tends overall to reduce summertime precipitation, especially over the west side of Central America (Wang 2007; Amador 1998; Mo et al. 2005; Cook and Vizy 2010).

To monitor the correlation between the CLLJ and SRI, we introduce, following the definition of (Wang 2007), the CLLJ index as the area-averaged MJJAS 925 hPa zonal wind in the region represented in Fig. 2a (12.5∘N–17.5∘N, 80∘W–70∘W) multiplied by − 1. This is justified by the fact that easterly winds feature a maximum (larger than 13 m/s) in the lower troposphere at about 925 hPa (Fig. 2a). The 5-year running mean of the CLLJ index is shown in Fig. 2b. The correlation between CLLJ and SRI in summer is − 0.58 in observations, with models generally reproducing such relationship well (Table 2).

Fig. 2
figure 2

Anomalously strong Caribbean low-level jet and tropical Atlantic-tropical Pacific temperature difference. a 2015–2019 zonal wind anomaly at 925hPa (dataset: ERA5). Contours show the climatological mean (1979–2019) of the zonal wind component, which highlights the position of the Caribbean low-level jet. b Five-year running mean of the Caribbean low-level jet index (Wang 2007) (dataset: ERA5). c 2015–2019 SST anomaly relative to the 1921–1970 climatological mean (contours). The long term positive trend is removed to better highlight the interbasin temperature difference (dataset: ERSSTv5). d Five-year running mean of the the Tropical Pacific-Tropical Atlantic SST difference index (Fuentes-Franco et al. 2015). Dataset: ERSSTv5

Full size image

Interannual variability of precipitation over Central America is affected by interactions with Pacific and Atlantic Ocean SSTs. Warm Atlantic-cool Pacific (cool Atlantic-warm Pacific) conditions cause a weakening (strengthening) of the CLLJ and favor increased (decreased) precipitation over Central America (Taylor et al. 2002, 2011; Fuentes-Franco et al. 2015, see also Fig. S6).

To quantitatively define the relationship between SRI and the Tropical North Atlantic (TNA)-Tropical North Pacific (TNP) SST difference, we calculate the TNATNP SST index (TNATNP hereafter) as defined by Fuentes-Franco et al. (2015). The TNATNP is the difference between the TNA SST averaged in the region between 5 to 22∘N and from 85 to 35∘W, and the TNP SST averaged in the region between 9 to 27∘N and from 110 to 90∘W (Fig. 2c). As expected, the TNATNP index is positively correlated to SRI (0.6, Table 2), with models generally featuring an even larger correlation (between 0.7 and 0.8). The 5-year running mean of the TNATNP anomaly index indicates the persistence of conditions characterized by TNP warmer than the TNA (Fig. 2d).

Hereafter, we refer to the 2015–2019 Central American SRI negative anomaly (Fig. 1d) as “prec_event_1519”. Averaging all four observation datasets, we obtain a 5-year anomaly of − 37 mm/month, which is the largest (in magnitude) measured over the observational record (i.e., 1921–2019). Similarly, for the CLLJ and the TNATNP, we refer to the mean 2015–2019 MJJAS CCLJ and 2015–2019 MJJAS TNATNP anomalies as the “cllj_event_1519” and ‘tnatnp_event_1519”. The values of these anomalies are 14.8 m/s and − 0.23 K, respectively (Fig. 2b and c). The CLLJ has been persistently intense during 2014–2019, resulting in record high values of the 5-year running mean values of the CLLJ index (Fig. 2b). Likewise, the mean 2015–2019 SSTs in the tropical Pacific were generally higher than the SSTs in the tropical Atlantic by about 1 K (Fig. 2c), resulting in large negative values of TNATNP (Fig. 2d).

2.5 Estimate of probability and risk ratios

The probability of occurrence of prec_event_1519 is estimated according to the empirical probability distribution of the 5-year SRI anomalies in the conterfactual climate, i.e., a climate undisturbed by anthropogenic influence (e.g., NAS 2016; Hauser et al. 2017). As noted by Hauser et al. (2017), choices of different counterfactual climates (e.g., pre-industrial vs. early twentieth century) can lead to contradicting results. We account (when possible, see Table 1) for three counterfactual climates: (1) CTRL, (2) NATURAL, and (3) PAST. PAST employs historical simulations forced with observed boundary conditions (i.e., HIST) over time period (1921–1970) when the anthropogenic impact on climate was substantially smaller than today.

As in Pascale et al. (2020), to build the CTRL empirical probability distribution, we randomly selected a 50-year and 5-year sequence (nonoverlapping) and then calculate the anomaly of the 5-year period relative to the 50-year climatology. This choice mimics the 2015–2019 mean minus the 1921–1970 mean. We then repeat this process 10,000 times to form a distribution of the 5-year SRI anomalies. We deal with the NATURAL runs in the same way after concatenating all the ensemble members. For the HIST runs, we concatenate the ensemble members after taking years 1921–1970.

To evaluate the decadal change in the probability of occurrence during the historical, present (factual climate) and future projected climate, we empirically estimate a decadal-varying probability distribution of the 5-year MJJAS rainfall anomalies (relative to the 1921–1970 mean) using the ALLFORC experiments (HIST and SSP5-8.5 or SSP2-4.5 in the case of SPEAR_MED). Following the method used in Pascale et al. (2020), we estimate the empirical probability distribution for a 20-year time window centered around the referenced year (so, for example, decade 2010 is built from all years from 2001 to 2020). This choice is a trade-off between a time period not too wide (in order to assume the stationarity of the probability distribution) and a number of instances large enough to allow for sufficiently accurate estimates of probabilities of rare events (e.g., 100-y return time). Once we have the decadal probability distribution, we can estimate the probability of occurrence, for each 20-year period, of prec_event_1519 for any random 5-year segment within the 20-year time window. We repeat this calculation every 5 years, that is, 1921–1940, 1926–1945, ..., 2076–2095, 2081–2100, to obtain a time-varying estimate of the probability.

We estimate the risk ratio of prec_event_1519 by computing the ratio of the probability that the MJJAS precipitation deficit is below the chosen threshold in the factual world and the probability that the MJJAS precipitation deficit is below the same threshold in the counterfactual world. For the probability of the factual world (at the time of the event), we take a weighed mean of the probabilities for the 2015 (i.e., 2006–2025) and 2020 (2011–2030) decades. The 95% confidence intervals in these probabilities were estimated by applying bootstrap-with-replacement resampling. The same approach described so far is applied to the cllj_event_1519 and tnatnp_event_1519 events.

3 Results

3.1 Recent 40-year trend

The recent negative rainfall trends are common to most of Central America (Fig. 3a and Fig. S2). Previous studies (Neelin et al. 2006) found that simulations from climate models of the CMIP3 archive (Meehl et al. 2007) generally show a good intermodel agreement on this pattern. This led Neelin et al. (2006) to wonder whether this regional precipitation anomaly is a consequence of ACC. After analyzing CMIP3 (Meehl et al. 2007) models preindustrial runs, they concluded that attribution to ACC of the observed negative trend is plausible but inconclusive.

Fig. 3
figure 3

Multi-decadal drying trend in Central America. a 1979–2019 trend in MJJAS precipitation totals (GPCC dataset; Schneider et al. 2013). b Probability distribution of 40-year SRI trends from unforced, long-term CTRL runs (see Table 1). The observed SRI negative trend, with its 95% significance confidence interval, is shown by the black vertical line and grey bar. c Ensemble distribution of the 1979–2019 SRI trends. The trend ensemble mean, for each ensemble, is shown by the color marks on the upper x-axis. The pale blue denotes the range of the forced (i.e., ensemble mean) trends. d Ensemble distribution of the 1979–2019 SRI trends from the ALLFORC and NATURAL (Table 1) ensembles. Precipitation trends are estimated using least-squares linear trends for each grid point. Significance of trends at the 5% level is assessed through a t-test

Full size image

Here, we first evaluate how likely the 1979–2019 SRI trend is in a counterfactual climate with no anthropogenic influence, i.e., in the long term control runs (CTRL, Table 1). The CTRL 40-year trend distributions (Fig. 3b) suggest that this is indeed a rare yet possible event in four out of five model ensembles analyzed in this study (0.9% in SPEAR_MED, 0.6% in SPEAR_LO, 0.07% in FLOR_FA, 0% in CESM and 1.1% in MPI-GE). Differences in these probability estimates might depend on the model’s ability to reproduce the variety of trends due to internal variability. For example, CESM largely underestimates natural variability of Central American summer rainfall compared to observations (standard deviation of 25 mm/month vs 33 mm/month). Nevertheless, we can state that a 40-year trend like that observed recently is likely to be possible even without any anthropogenic influence on the climate.

As a second step of our analysis, we then analyze the 1979–2019 trends from the ALLFORC ensembles (Table 1), whose forcing includes both natural and anthropogenic effects. Fig. 3c shows the distribution of SRI 1979–2019 trends for each model ensemble. These distributions are centered around their forced value, with their widths determined by internal climate variability. For each model, the trend’s ensemble mean, which represents the forced signal, is close to zero (between − 0.3 and 0 mm/month/year) and substantially smaller than observations (− 1.3 mm/month/year). This indicates a small forced signal relative to the noise generated by internal variability. Large negative trends comparable to the observed trend are achieved only by a few ensemble members (Fig. 3c).

Finally, for models for which we have the NATURAL experiments (SPEAR_MED and FLOR_FA), we compare the trend distributions across the ensemble members for the ALLFORC and NATURAL experiments (Fig. 3d). This can be very informative if we recall that ALLFORC and NATURAL differ only by the presence of anthropogenic forcing. Thus, any significant difference in the two distributions is to be attributed to anthropogenic forcing. We find that, for both SPEAR_MED and FLOR_FA, the two distributions are qualitatively very similar (Fig. 3d). A Kolmogorov-Smirnov test performed on the two arrays of trends, for each model, reveals that the probability that the two arrays are drawn from the same distribution is 94% for both models. This test so does not pass the 5% level, although it falls close to it (i.e., 95%). This results so suggests, at least for these two models for which we have NATURAL, that the addition of anthropogenic forcing to natural forcing does not lead to a substantial difference in the 40-year trends distribution across the ensemble.

The overall evidence presented so far suggests that the observed recent drying trend may be mostly due to internal, multi-decadal variability rather than ACC. However, we admit that we cannot totally exclude the role of ACC. This is because we are unable to make a more precise statement (of the form “ACC has increased the likelihood of the 1979-2019 SRI trend by X times”) given the scarcity of points, thirty, used to reconstruct the probability distribution of the 1979–2019 SRI trends (Fig. 3d). Instead, we would have needed hundreds to thousands of equivalent (i.e., under the same 1979-2019 natural only and natural plus anthropogenic forcing) estimates of the 1979–2019 SRI trends. These would have then allowed us to reconstruct reliably the probability distribution of the 1979–2019 SRI trends and so estimate the probability of having a 40-year trend smaller (i.e., more negative) than the one observed in 1979–2019 in both the factual and counterfactual climate.

3.2 Attribution of the 2015–2019 drought event

Figure 4 shows, for each model ensemble, the risk ratio of prec_event_1519. These values are obtained by the ratio of the probabilities of occurrence of prec_event_1519 in the factual world and the probability of occurrence of prec_event_1519 in three different counterfactual climates (CTRL, NATURAL and PAST, see Section 2.5). For SPEAR_MED, we have two different sets of values, depending on whether years 2015–2030 are drawn from SSP2-4.5 or SSP5-8.5 (Fig. 5a).

Fig. 4
figure 4

Risk ratios for prec_event_1519 in 2015–2019 (PRES) for all ensembles and for all the available counterfactual climates (PAST, CTRL, NATURAL). Bars show the median and the 95% uncertainty intervals on a logarithmic axis. NATURAL runs are not available for the CESM and MPI-GE. A simple average is used to synthesize the results

Full size image
Fig. 5
figure 5

Probability of occurrence of a 5-year MJJAS rainfall anomaly equal to or worse than in 2015–2019 in HIST, SSP2-4.5, and SSP5-8.5 in SPEAR_MED. Shading denotes the 95% confidence interval from bootstrap resampling. The red constant line denotes the CTRL probability for such an event, and the red constant dashed line is from the NATURAL run after concatenating all 30 ensemble members

Full size image

All ensembles, with the exception of SPEAR_MED ALLFORC4.5, indicate a risk ratio significantly larger than one for all three counterfactual climates. The difference between SPEAR_MED ALLFORC4.5 risk ratio and SPEAR_MED ALLFORC8.5 risk ratio originates from years 2015–2030, which lead to significantly different values of the probability of occurrence of prec_event_1519 (Fig. 5a). These differences may be in part explained by the differences in radiative forcings between SSP2-4.5 and SSP5-8.5, which are minimal until 2020–2030 and start differing substantially only afterwards (Riahi et al. 2017). Differences may be also attributable to the limited number of ensemble members (30), which does not allow a complete quantification of internal variability (Fig. 5a).

For the second part of the twenty-first century, the probabilities of occurrence of prec_event_1519 in SSP2-4.5 and SSP5-8.5 diverge substantially (Fig. 5b). In SSP2-4.5 (middle of the road pathway) it remains below 20% by the end of this century. Given a probability of around 2% of prec_event_1519 in the counterfactual climates, this implies a risk ratio which remains below 10. In SSP5-8.5 (high emission scenario), this probabilities increases swiftly, reaching 90% by the end of the century, that is a risk ratio of nearly 40. SPEAR_LO and FLOR_FA share similar values (median values between 2 and 5). For MPI (100 ensemble members), we use RCP8.5 for years 2005–2030, but we find no significant differences when RCP4.5 and RPC2.6 are instead used (not shown).

Combining all different risk ratios for different models, using a simple average to synthesize the results, the probability of an event like the observed 2015–2019 meteorological drought has increased by a factor of 4 relative to the two counterfactual climates PAST and CTRL, with confidence intervals (3, 7) and (3.5, 5) respectively, and by a factor 1.5 (0.95, 2) relative to NATURAL. Let us note, however, that the value for NATURAL is obtained averaging only three cases, so we should bear in mind that is a less comprehensive and somewhat biased result. Overall, results summarized in Fig. 4 suggest that ACC has significantly increased the likelihood of a multi-year drought like that occurred in Central America between 2015 and 2019.

3.3 The role of the CLLJ and Tropical Atlantic-Tropical Pacific temperature SST differences.

The 2015–2019 meteorological drought unfolded in association with a very strong CLLJ (Fig. 2a) and a negative TNATNP (Fig. 2c). In particular, the 5-year mean of the CLLJ index features the largest positive anomalies over the entire ERA5 time span (Fig. 2b) while the 5-year mean of TNATNP features a negative value equalled only other two times in the last hundred years (Fig. 2d).

The observed mean 2015–2019 925 hPa zonal wind anomalies and SST anomalies have a strong resemblance to those projected for the end of this century under the highest emission scenario (Fig. 6). Strengthening of the CLLJ and faster warming up of the tropical Pacific than the tropical north Atlantic is indeed a robust feature across climate models (Rauscher et al. 2008, 2011; Fuentes-Franco et al. 2015). It is therefore natural to ask whether the likelihood of the unusual 2015-2019 CLLJ and TNATNP anomaly was increased by ACC. To address this question, we estimate the probability of occurrence of cllj_event_1519 and tnatnp_event_1519 (Section 2.4) in the factual and counterfactual climates. We perform this analysis only for SPEAR_MED and SPEAR_LO, for which we have data readily available. However, given the similar future projections in terms of CLLJ and SSTs, we speculate that similar results would be obtained from other ensembles too. The CLLJ index features a long term significant trend in both NATURAL and HIST (SPEAR_MED). Since this trend is not seen in reanalyses, and there is no long-term forcing in NATURAL which can physically justify it, it must be a model drift unrelated to anthropogenic forcing. We therefore subtract from the CLLJ index the trend estimated from NATURAL to both NATURAL and HIST before estimating decadal probabilities of occurrence and risk ratios.

Fig. 6
figure 6

SPEAR_MED projected strengthening of the Caribbean low-level jet and tropical Atlantic-tropical Pacific temperature difference. a Projected (2071–2100 vs. 1921–1970) June-July anomalies in the zonal component of the 925 hPa zonal wind (black contours show the climatological easterly component) under SSP5-8.5. b Projected (2071–2100 vs. 1921–1970) MJJAS SST anomalies (black contours show the climatological mean) under SSP5-8.5

Full size image

Figure 7 shows the risk ratio of cllj_event_1519. We have contrasting results across different counterfactual climates (e.g., in SPEAR_LO) and for the same model (i.e., SSP2-4.5 vs. SSP5-8.5). Overall, there is not strong evidence for the 2015–2019 CLLJ mean anomaly to be attributable to ACC. In SPEAR_MED (SSP5-8.5), risk ratios of cllj_event_1519 are close to one (indicating no attribution to ACC) even though those of prec_event_1519 are clearly larger than one (Fig. 4). Similarly, in SPEAR_LO, risk ratios for prec_event_1519 are clearly larger than one, but risk ratios for cllj_event_1519 are larger than one only when PAST is chosen as the counterfactual climate. This means that for almost all cases in which prec_event_1519 is attributable to ACC, the 2015-2019 CLLJ anomaly is not attributable to ACC. While the CLLJ has not shown fluctuations unambiguously attributable to ACC during 2015–2019, the decadal evolution of the probability of cllj_event_1519 in SPEAR_MED suggests that this will be the case in the next few decades (Fig. 9a and b).

Fig. 7
figure 7

Risk ratios for cllj_event_1519 in 2015–2019 (PRES) for all SPEAR ensembles and for all the available counterfactual climates (PAST, CTRL, NATURAL). Bars show the median and the 95% uncertainty intervals on a logarithmic axis. Risk ratio values for SPEAR_MED SSP2-4.5 and SSP5-8.5 differ because of the different probability for the event15_19 during the years 2015–2030 used in the 20-year time window

Full size image

Figure 8 shows the risk ratios associated with tnatnp_event_1519. Results here are much more consistent across the SPEAR models and with those about the prec_event_1519 event, and they indicate that ACC has increased the likelihood of the tnatnp_event_1519 by a factor 2 (with a 95% confidence interval between 1.5 and 5). The NATURAL data are lacking for this case, although those do not generally lead to very differ results from those of, e.g., CTRL (Fig. 7) and therefore that does not alter the validity of our conclusions.

Fig. 8
figure 8

Risk ratios for tnatnp_event_1519 in 2015–2019 (PRES) for all SPEAR ensembles and for all the available counterfactual climates (PAST, CTRL). Bars show the median and the 95% uncertainty intervals on a logarithmic axis

Full size image

The emergence of the ACC impact on TNATNP is faster than that on CLJJ, as evident from the time evolution of the probability of occurrence (e.g., Fig. 9a vs. c, b vs. d). We conjecture that may be due to the larger level of background noise of a dynamical, atmospheric field as the CLLJ as compared the more thermodynamically constrained SST field. This leads us to speculate that—differently from SST, for which ACC is already detectable (Chan and Wu 2015)—the signal associated with ACC on the regional circulation, specifically the CLLJ, has not yet emerged from the noise of internal variability. For example, previous work on the strengthening of the North Atlantic Subtropical High—which is connected to the CLLJ, although it is not the only driver of the CLLJ—is also inconclusive about the time of emergence of the ACC forced signal (Li et al. 2011, 2012, 2013, Diem 2013).

Fig. 9
figure 9

a–b Probability of occurrence of a 5-yr MJJAS CLLJ index anomaly equal to or worse than in 2015-2019 in HIST, SSP2-4.5, and SSP5-8.5 in SPEAR_MED. Shading denotes the 95% confidence interval from bootstrap resampling. The red constant line denotes the CTRL probability for such an event, and the red constant dashed line is from the NATURAL run after concatenating all 30 ensemble members. c–d As for a–b but for the TNATNP index

Full size image

4 Discussion and conclusions

In this study, we have analyzed the recent drying trend observed from the late 1970s and the 2015–2019 Central American meteorological drought to determine whether these two events have been made more likely by ACC.

From the comparison of large ensembles driven by natural forcing only and natural plus anthropogenic forcing, we conclude that the observed 1979–2019 drying trend is consistent with internal, multi-decadal variability, although our analysis on this was not conclusive and therefore we cannot totally exclude the role of ACC. Instead, we found strong evidence that the 2015–2019 rainfall deficit was made four times more likely by ACC (Fig. 4). Closely related to that, we also found that ACC increased the probability (1.5 to 3), for the period 2015–2019, of large, negative SST differences between the Tropical North Atlantic and Tropical North Pacific, which are known to drive meteorological droughts in Central America.

At this point it is worth speculating on why we managed to attribute the 2015–2019 rainfall deficit to ACC while we did not for the 1979–2019 trend. First, as already discussed in Secttion 3.1, we are unable to quantify a risk ratio for the 1979–2019 trend event—and thus quantify by how many times ACC made it more or less likely—because we cannot empirically reconstruct the probability distribution for this event due to scarcity of data. Second, we note that the effect of ACC starts becoming noticeable on the 5-year SRI distribution’s tail only after 2010 in SPEAR_MED (Fig. 5a) and after 2000 in FLOR_FA. Therefore, while this allows us to see the effects of this climatic change on a short-term event like prec_event_1519, it may not directly lead to a detectable effect on the 1979–2019 trend, whose estimate is based, for the most part, on years during which the effect of ACC was not yet detectable (at least according to these two models).

The extent of the precipitation anomaly observed between 2015 and 2019 revealed a key role for rainfall in this multi-year drought, leading us to focus on precipitation as the dominant variable to characterize this event. Indeed, recent studies highlight the projected precipitation reduction over Central America as the key factor in the increase of drought risk in the coming decades (Giorgi 2006; Rauscher et al. 2008; Fuentes-Franco et al. 2015; Almazroui et al. 2021; Depsky and Pons 2020). However, it is also important to mention that regional ongoing warming trends (Aguilar et al. 2005; Gourdji et al. 2015) play an important role in driving enhanced evapotranspiration and thus further exacerbating future droughts (Hidalgo et al. 2021, 2019; Alfaro-Cordoba et al. 2020). An example of that is the 2013–2016 Pan-Caribbean drought, which was not driven by a precipitation deficit but it was a direct consequence of exceptionally high evapotranspiration due to higher temperatures (Herrera et al. 2018). Furthermore, higher surface temperature from surrounding oceans along with drier land conditions can induce stronger lower-to-middle tropospheric stability, especially in the early summer, by increasing the midtropospheric dry static energy relative to near-surface moist static energy (Seth et al. 2011, 2013; Giannini 2010). Increased lower-to-middle tropospheric stability suppresses convection, reducing precipitation in monsoonal regimes where summertime precipitation is mostly associated with convective systems (Pascale et al. 2017, 2019). Higher temperatures therefore may exacerbate the consequences of a meteorological drought in multiple ways. Future work should include analysis of the risk associated with compound events in Central America driven by exceptionally high temperatures and prolonged rainfall deficits.

Although a multi-year meteorological drought like the one that hit Central America between 2015 and 2019 is a rare event, anthropogenic climate change has significantly increased the likelihood for such a prolonged rainfall deficit to occur. All of the models analyzed in this study suggest that the probability of occurrence of this event will continue to increase into the future, with a rate depending on the specific socioeconomic scenario (Eyring et al. 2016; O’Neill et al. 2017; Riahi et al. 2017). As a consequence, we expect that prolonged multi-year periods of below-average precipitation in the region will be more severe in the coming decades.

The more frequent occurrence of multi-year meteorological droughts in Central America shown in this study, along with increased aridity due to higher temperatures and the likelihood of increased water demand due to a growing population (Hidalgo 2021), suggests that managing water resources will likely be increasingly challenging in the decades to come. Advanced planning for drought years and measures to increase water efficiency for agricultural and urban use is therefore urgently needed at the national and regional level to better adapt to the new hydroclimatic conditions.